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HIGHLIGHTS 


•  An  AA-CG-MD  method  is  developed  to  simulate  the  fabrication  process  of  the  porous  anode  of  SOFC. 

•  This  method  is  employed  to  reconstruct  the  microstructure  of  the  porous  anode  vividly. 

•  Relevant  thermal  properties  are  calculated  in  a  rapid  and  accurate  manner  on  the  basis  of  the  method. 

•  Analyses  and  predictions  are  proceeded  to  capture  the  good  performance  for  the  porous  anode  of  SOFC. 
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Incorporation  of  nanoscale  catalysts  into  porous  structures  of  SOFC  has  been  proven  highly  successful  in 
increasing  active  sites  and  catalyst  utilization.  In  addition,  electrochemical  reactions  as  well  as  heat 
transfer  process  in  porous  anode  are  strongly  affected  by  complex  porous  structures.  It  is  believed  that 
study  of  anode  thermal  properties  are  critical  for  SOFC  design  and  operation.  In  this  work,  an  AA  model  is 
developed  for  nickel  and  YSZ  components  via  ASE,  and  a  CG  technique  is  further  applied  to  represent  Ni 
and  YSZ  beads  by  VMD,  which  are  then  self-assembled  to  capture  the  anode  nanostructure  via  GRO- 
MACS.  LAMMPS  is  then  employed  to  evaluate  average  thermal  properties  of  the  porous  anode.  It  is  found 
that,  at  low  Ni  content  (<30  vol%),  thermal  conductivity  increases  with  increasing  temperature  due  to 
lattice  vibrations.  Instead,  the  anode  exhibits  metallic  behavior  due  to  rich  nickel  phase.  Thermal 
expansion  of  the  anode  increases  with  increasing  nickel  content.  Average  thermal  properties  of  the 
anode  are  validated  by  open  literature  data  with  good  agreement.  This  approach  is  considered  to  be 
applied  to  analyze  nanostructures,  heat  transfer  and  temperature  distribution  in  the  porous  anode,  and  is 
also  useful  to  capture  thermal  performance  of  SOFC  and  stack. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cell  (SOFC)  is  an  attracting  energy-converting 
device  due  to  its  high  energy  efficiency  [1,2]  and  fuel  flexibility 
[3,4],  as  well  as  low  pollutant  emissions  [5].  SOFC  employs  multi¬ 
layered  electrodes  which  comprise  ceramic  and  metallic  materials 
with  different  thermal  properties.  All  components,  e.g.,  anode, 
cathode,  electrolyte,  etc.,  have  to  provide  a  well-adjusted  heat 
diffusion,  thermal  expansion,  mechanical  strength  and  so  on,  both 
at  material  interfaces  and  inside  materials.  The  performance  of 
SOFC  is  not  only  determined  by  intrinsic  material  properties,  but 
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also  subject  to  the  fabrication  technology  of  functional  layers. 
Basically,  that  can  be  improved  by  the  application  of  modified 
materials  using  appropriate  technologies.  It  is  well  known  that  the 
higher  operating  temperatures  in  SOFC  improve  its  performance  in 
comparison  with  low-temperature  fuel  cells,  e.g.,  PEMFC,  DMFC, 
PAFC,  etc.  In  addition,  high  operating  temperatures  make  it  feasible 
to  carry  out  electrochemical  reactions  in  the  active  sites  of  the 
porous  anode  with  cheap  catalysts.  On  the  other  hand,  the  internal 
reforming  reactions  of  hydrocarbon  fuels,  e.g.,  methane,  are 
contributing  to  keep  a  good  balance  of  heat  generation  and  con¬ 
sumption.  Heat  is  mainly  produced  in  the  exothermically  electro¬ 
chemical  reactions.  However,  it  is  noticed  that  too  much  cooling 
caused  by  the  endothermic  steam  reforming  reactions  has  a 
negative  effect  on  the  output  performance  of  SOFC  [6],  which  is 
remarkably  opposed  to  desired  thermal  management  in  the  porous 
anode  for  SOFC.  Thus,  studying  heat  transfer  in  order  to  effectively 
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enhance  the  thermal  management  inside  the  porous  anode  has 
been  attracting  significant  attention  in  the  recent  years.  This  means 
that  to  correctly  predict  the  thermal  properties  of  the  porous  anode 
of  SOFC  is  important  in  developing  high-performance  SOFCs. 

In  order  to  generate  high  power  density  and  achieve  reaction- 
related  efficiency  as  high  as  possible,  abundant  triple-phase 
boundaries  (TPBs)  [7,8],  where  the  electrochemical  reactions  take 
place,  should  be  available  in  the  porous  electrodes.  This  explicitly 
represents  a  principle  that  to  sinter  high  porosity  in  the  functional 
anode  seems  extremely  important  for  the  extensive  distribution  of 
so-called  TPBs.  Unfortunately,  for  the  well-dispersed  porous  anode 
from  nanoscale  particles,  e.g.,  Ni  and  yttria-stabilized  zirconia  (YSZ) 
[9,10],  etc.,  significant  reduction  of  the  thermal  conductivity  is  ex¬ 
pected  due  to  discontinuities  inside  mediated  materials,  which  is 
described  adequately  by  effective  medium  theory  (EMT)  [11].  The 
mechanical  strength  of  the  porous  anode  is  also  prone  to  be  low¬ 
ered  to  a  certain  extent.  It  is  thus  extremely  important  to  keep  a 
good  balance  among  high  porosity,  considerable  thermal  diffusion 
and  acceptable  mechanical  strength  in  the  porous  anode  of  SOFC.  In 
addition,  high  temperature  gradients  will  easily  take  place  along 
and  normal  to  the  direction  of  fuel  flow  due  to  mass  diffusion  or 
concentration  drop  of  the  fuel  gas.  This  will  be  mainly  attributed  to 
the  unreasonable  cooling  effect  due  to  the  high  supply  rate  of  the 
air  flow  in  the  cathode.  It  is  found  that  a  sharp  temperature 
gradient  or  thermal  stress  is  primarily  responsible  for  cracking  in 
the  porous  anode,  even  in  the  electrolyte.  The  temperature  gradient 
or  the  thermal  stress  accounted  for  by  heat  transfer  inside  the  SOFC 
can  be  computed  numerically  but  really  depends  on  specific  heat 
capacity  as  well  as  heat  flow  rate  and  thermal  conductivity  (TC),  etc. 

In  this  work,  an  all-atom  coarse-grained  molecular  dynamics 
(also  named  AA-CG-MD)  method  is  employed.  In  brief,  a  target 
domain  with  cubic  50  nm  (typically  limited  by  GROMACS,  the  code 
used  to  reconstruct  the  nanostructures  of  the  porous  anode, 
especially  the  TPB  region)  will  be  implemented  from  the  ab  initio 
stage,  e.g.,  atoms,  molecules,  lattices,  etc.  The  box  will  be  filled  with 
beads,  also  named  “pseudo  atoms”,  representing  groups  of  atoms/ 
molecules  to  be  computed  in  the  MD  approach.  The  model  may  be 
developed  for  the  reforming  and  electrochemical  reactions,  mass 
and  charge  transfer,  and  heat  diffusion.  The  predicted  material 
properties  in  this  study  based  on  the  open  literature  data  [12-16] 
are  linked  to  the  compositions  and  nanostructures  of  the  mate¬ 
rials.  Some  basic  principles  govern  the  choices  of  SOFC  materials. 
For  example,  the  electrodes  should  have  high  electrical  conduc¬ 
tivity  for  the  charge  transport,  high  catalytic  activity  for  the 
involved  reforming  and  electrochemical  reactions,  adequate 
porosity  for  the  gas/vapor  diffusion,  good  physical/chemical 
compatibility  with  the  electrolyte  and  interconnect  for  the  long¬ 
term  stability. 

As  experimental  measurements  tend  to  be  more  costly  and 
difficult  to  carry  out  compared  with  the  simulation  work,  which  is 
speedy  and  relatively  easy  once  the  model  is  programmed  and 
validated,  their  ability  to  investigate  the  effects  of  different  material 
parameters  and  operating  conditions  on  the  performance  of  SOFC  is 
critically  limited.  Numerical  modeling  and  simulations  incorpo¬ 
rating  the  known  physical  and  chemical  phenomena  occurring  in¬ 
side  SOFC  to  predict  its  performance  are  quite  important  for  the 
understanding  and  technological  improvement  of  SOFC.  Though 
there  is  no  shortage  of  research  on  the  numerical  modeling  for 
SOFC,  the  analysis  focusing  on  the  effects  of  nanostructures  in  the 
porous  anode  on  thermal  properties  is  relatively  rare  in  the  open 
literature.  Therefore,  it  seems  that  such  a  research  effort  is  quite 
essential  for  improving  the  performance  of  SOFC.  In  this  case,  the 
thermal  properties,  electrochemical  performance  and  their  inter¬ 
play  are  supposed  to  be  taken  into  account  to  properly  extract  the 
best  performance  for  SOFC  to  practical  applications.  This  work 


describes  a  detailed  nanoscale  model  for  the  numerical  simulations 
in  the  porous  anode  of  SOFC  with  a  representative  Ni/YSZ-YSZ- 
LSM/YSZ  setup  [17,18].  The  thermal  properties  of  the  porous  anode 
are  determined  and  evaluated  against  the  open  literature  data  [19- 
22].  Anodic  nanostructure,  as  well  as  effects  of  Ni  contents,  and 
sintering  conditions  of  the  anode  on  its  thermal  properties  are 
systematically  investigated.  Finally,  the  desired  material  composi¬ 
tions  and  anodic  nanostructures  are  deduced. 

2.  Modeling  development 

In  this  study,  the  AA-CG-MD  approach  is  firstly  implemented  to 
reconstruct  microscopic  structures  for  the  porous  anodes  at  the 
nanoscale,  and  the  obtained  structures  are  further  applied  to 
calculate  thermal  properties  of  the  porous  anodes,  e.g.,  average  TC, 
thermal  expansion  coefficient  (TEC),  volumetric  heat  capacity 
(VITC),  etc.  The  size  effect  on  TC  is  also  evaluated.  These  contribute 
to  further  predict  the  probability  of  thermal  failure  and  evaluate  the 
overall  performance  of  the  porous  anode  of  SOFC  [23-25]  by  the 
MD  modeling  [26].  The  precise  details  of  the  method  are  presented 
in  the  following  three  steps. 

2  A.  All-atom  modeling 

In  classical  MD  modeling,  a  single  potential  energy  surface  is 
represented  in  the  force-field,  which  is  a  consequence  of  Born— 
Oppenheimer  (BO)  approximation  [27,28].  Therefore,  it  is  possible 
to  compute  the  wave-function  and  the  energy  of  an  atom/molecule 
in  finite  and  less  complicated  steps. 

^total  =  ^electronic  x  0 nuclear  O) 

where  W  referring  to  Equation  (1),  allows  the  wave-function  of  an 
atom/molecule  to  be  broken  into  its  electronic  and  nuclear 
components. 

In  excited  states,  when  electrochemical  reactions  or  a  more 
accurate  representation  are  needed,  electronic  behavior  can  be 
obtained  from  first  principles  [29]  by  using  a  quantum  mechanical 
method,  such  as  Density  Functional  Theory  (DFT)  [9,20].  This  is 
known  as  the  AA  modeling  [30].  Due  to  treating  electronic  degrees 
of  freedom,  the  computational  cost  of  this  method  is  much  higher 
than  classical  MD,  which  implies  that  the  AA  modeling  is  limited  to 
smaller  systems  and  shorter  periods  of  time.  In  general,  this 
modeling  is  usually  made  in  the  close  neighborhood  of  the  atom / 
molecule-based  systems.  Although  various  approximations  may 
be  used,  these  are  based  on  theoretical  considerations  instead  of 
empirical  fitting.  The  AA  calculations  generate  a  huge  amount  of 
information  that  is  not  available  from  empirical  methods,  such  as 
density  of  electronic  states  or  other  electronic  properties.  A  sig¬ 
nificant  advantage  of  using  the  AA  modeling  is  the  ability  to  study 
nanostructures,  even  interactions  that  involve  breaking  or  forma¬ 
tion  of  covalent  bonds,  which  correspond  to  multiple  electronic 
states  and  potential  energy.  In  this  study,  the  AA  modeling  is  mainly 
applied  to  generate  the  nanocrystals  for  sintering  powders  at  the 
atom-scale,  i.e.,  Ni  and  YSZ  particles.  The  face-centered  cubic  (FCC) 
nickel  nanocrystals  are  generated  according  to  the  AA  modeling, 
and  partially  stabilized  tetragonal  YSZ  nanocrystals  are  also  pro¬ 
duced  very  well  by  the  same  technique,  and  then  both  of  them 
compose  the  so-called  “green  body”  (GB)  [31]  for  sintering  the 
porous  anodes.  Although  no  open  literature  data  supports  that  the 
crystalline  nanostructure  shows  higher  oxygen-ion  conductivity 
than  the  amorphous  nanostructure  in  the  porous  anode,  extensive 
research  work  has  been  focused  on  this  over  the  past  decades.  It  has 
been  found  that  the  crystalline  nanostructure  is  more  stable  ener¬ 
getically  compared  to  the  amorphous  nanostructure  in  the  porous 
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anode  [32].  For  these  smaller  systems,  the  AA  modeling  re¬ 
constructs  the  crystalline  nanostructures  of  the  sintering  particles 
for  the  porous  anode  in  excellent  agreement  with  SEM/TEM  [22,33] 
patterns,  and  to  some  extent,  for  much  larger  systems,  this 
approach  can  still  reveal  the  approximate  morphological  occu¬ 
pancy  for  Ni/YSZ  composition,  although  the  interactions  of  Ni/YSZ 
particles  become  poorly  determined. 


2.2.  Coarse-graining  modeling 

On  the  other  end  of  nanoscale,  there  are  coarse-grained  and 
lattice  models.  Instead  of  explicitly  representing  every  atom/ 
molecule  of  the  system,  one  uses  “pseudo-atoms”  to  represent 
groups  of  atoms.  Further  speaking,  CG  [34,35]  is  composed  of 
replacing  an  atomic  description  of  an  atom/molecule  with  a  lower- 
resolution  coarse-grained  particle  that  averages  or  smoothes  away 
fine  details. 

The  MD  simulations  on  very  large  systems  may  require  such 
huge  computing  resources  that  they  cannot  easily  be  studied  by  the 
traditional  AA  modeling.  Similarly,  simulations  at  long  time-scales 
(typically  more  than  1  ps)  are  prohibitively  expensive,  because 
they  require  so  many  time-steps.  In  these  cases,  one  can  sometimes 
handle  the  problem  by  using  reduced  or  simplified  representations, 
also  named  the  coarse-grained  (CG)  models,  which  have  been 
developed  to  investigate  many  complex  systems,  e.g.,  lipid  mem¬ 
branes,  proteins,  ceramics,  etc.,  at  the  longer  time/length-scales. 
The  following  formula  calculates  the  system’s  potential  energy  as 
a  sum  of  individual  energy  terms  when  the  CG  representation  is 
implemented. 

^(rij)  =  ^bonded  (rij )+  ^nonbonded  (rij )  (2) 

where  the  components  of  the  bonded/nonbonded  contributions 
are  specified  by  the  following  equations,  as  outlined  separately  in 
Equations  (3)  and  (7). 

^bonded  (rij)  m  ^bond  (Gj)  +  ^angle  (^ijk)  +  ^dihedral  (^ijkl)  (3) 

in  which, 

Wnj)  =  E^u(riJ-bu)2  (4) 

pairs 

Wangle  (^ijk)  =  2^k  (^ijk  _  ^ijk)  ^ 

angles 

^dihedral  (Vijkll  =  ^ijkl  +  C0S  f^ijkl  —  ^ijkl)  j  (®) 

dihedrals 

where  ry  is  the  distance  between  bead  i  and  bead  j;  0yk  is  the  angle 
among  beads  i,  j  and  k;  <pyki  is  the  proper  dihedral  among  surface  ijk 
and  surface  jkl;  I<  represents  the  force  constant;  b°,  6°  and  <p°  are 
critical  values.  In  this  work,  the  bond  and  angle  terms  are  modeled 
as  harmonic  potentials  centered  equilibrium  bond-length  values. 
Empirically,  the  force  constants  of  JCbond,  Wangle  are  specified  as 
1250  kj  mol-1  nm-2,  25  kj  mol-1  rad-2  35],  respectively.  The 
involved  critical  values  are  determined  by  shape-based  coarse- 
graining  (SBCG)  tool  of  Visual  Molecular  Dynamics  (VMD)  package 
when  the  nanostructure  transforming  is  implemented  from  the  AA 
representation  to  the  CG  beads. 

In  contrast  with  the  bonded  terms,  it  is  much  more  computa¬ 
tionally  costly  to  calculate  the  nonbonded  terms  because  a  given 
bead  is  bonded  to  only  a  few  of  its  neighbors,  but  interacts  with 
every  other  bead  in  the  target  system. 


^nonbonded  (rij )  —  ^vanderWaals  (rij )  +  K 


electrostatic 


(7) 


As  one  of  the  nonbonded  interactions,  the  van  der  Waals  term 
falls  off  rapidly  with  distance.  It  is  typically  modeled  as  the 
Lennard-Jones  form,  as  treated  also  in  this  study. 


Wj(ru)  =  5>U 


(8) 


where  ey  is  the  depth  of  the  potential  well;  cry  is  the  finite  distance 
at  which  the  inter-particle  potential  is  zero,  and  a  cut-off  radius  of 
0.43  nm  [35]  is  employed  simultaneously  in  this  study. 

In  addition,  the  Coulomb  form  is  used  to  address  the  electro¬ 
static  term  of  the  nonbonded  interactions  in  Equation  (7). 

^Coulomb  (rij)  =  (9) 

where  q\  and  q j  are  the  charges  of  the  particles  i  and  j,  respectively; 
so  is  the  electric  permittivity  (8.85E-12  F  m-1). 

It  is  remarkably  worth  noting  that  those  coarse-grained  systems 
are  in  fact  limited  by  the  dynamical  accuracy  and  structure-based 
properties.  The  research  work  relevant  to  CG  is  in  its  infancy,  and 
its  theoretical  fundamental  basis  is  still  poorly  understood.  EIow- 
ever,  CG  has  already  been  widely  employed  in  the  area  of  MD. 

In  this  section,  as  sintering  powders,  the  nanocrystals  of  Ni  and 
YSZ  generated  by  the  above  AA  modeling  are  represented  by  a  few 
“pseudo-atoms”  which  gather  around  50  Ni  or  YSZ  atoms/mole¬ 
cules  into  each  “pseudo-atom”.  Subsequently,  the  virtual  bonded/ 
non-bonded  interactions  taking  place  among  the  “pseudo-atoms”, 
e.g.,  bonds,  angles,  torsion,  van  der  Waals  interactions,  etc.,  are 
extracted  as  well  via  the  CG  builder  (one  module  of  VMD  package), 
by  which  the  CG-MD  force-field  is  in  principle  achieved.  Unfortu¬ 
nately,  the  parameterized  process  for  CG  must  be  conducted  really 
dependent  on  experience,  by  matching  the  behavior  of  the 
modeling  to  appropriate  open-source  data  obtained  by  the  exper¬ 
iments  or  numerical  computations.  Ideally,  these  parameters 
should  account  for  both  bonded  and  non-bonded  contributions  to 
potential  energy  in  an  explicit  way.  Sometimes,  when  CG  is  done  at 
higher  levels,  the  accuracy  of  the  dynamic  description  may  be  less 
reliable  relatively.  But  a  few  coarse-grained  models  have  been  used 
successfully  to  tackle  a  wide  range  of  problems  in  structure-based 
research  work. 


2.3.  AA-CG-MD  modeling 


In  general,  MD  is  a  numerical  simulation  method  for  the  phys¬ 
ical/chemical  phenomena  in  complex  systems,  which  may  involve 
thermal  management,  mass  and  charge  transport,  together  with 
internal  reforming  reactions  and  electrochemical  reactions  in  the 
porous  anode  of  SOFC.  The  particles  (i.e.,  atoms,  molecules,  even 
“pseudo  atoms”)  from  the  large  systems  are  allowed  to  interact 
within  a  period  of  time,  and  then  give  a  view  of  the  motion  of  the 
computed  particles.  Theoretically,  these  are  determined  by 
numerically  solving  the  Newton’s  equation  of  motion  for  a  system 
of  interacting  particles,  as  below. 


Fj  =  m  j 


dj7j_ 

dt2 


aiT 


(10) 


where  m\  is  the  mass  of  the  particle  i,  r{  is  the  position  of  the 
particle,  V(ry)  is  the  potential  energy. 

The  forces  among  the  particles  and  the  potential  energy  (typi¬ 
cally  including  bond,  angle,  dihedral,  improper  dihedral,  van  der 
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Waals  and  electrostatic  potentials),  are  determined  by  the  prebuilt 
MD  force-fields,  e.g.,  MARTINI  [36]  (a  coarse-grained  force-field 
developed  by  Marrink  and  coworkers  at  the  University  of  Gronin¬ 
gen),  GROMACS  [37,38]  (the  force-field  optimized  for  GROMACS), 
VAMM  (Virtual  atom  molecular  mechanics),  etc.  Basically,  it  is 
almost  impossible  to  obtain  all  properties  of  such  complex  systems 
in  an  analytical  way  due  to  the  large  amount  of  particles  involved. 
Fortunately,  the  MD  simulation  appears  to  circumvent  perfectly 
this  limitation  by  utilizing  numerical  approaches,  if  it  is  ignored 
that  to  run  long  MD  simulations  appears  ill-conditioned  in  math¬ 
ematics.  This  phenomenon  is  mainly  attributed  to  cumulative  er¬ 
rors  during  integrations,  which  can  be  inhibited  slightly  with 
suitable  parameters  and  algorithms  implemented,  but  not  elimi¬ 
nated  entirely.  Despite  the  theoretical  support  behind  the  MD 
modeling  is  still  insufficient,  extensive  research  work  has  been 
focused  on  achieving  macroscopic  thermodynamical  properties  of 
the  target  system,  e.g.,  TC,  VHC,  TEC,  etc.,  in  terms  of  ergodic  hy¬ 
pothesis  (EH)  [39,40]  that  the  statistical  ensemble  average  is 
equivalent  to  the  time  average  on  the  simulated  system. 

In  this  part,  on  the  basis  of  the  above  AA  modeling  and  CG 
representations,  not  only  are  nanostructures  of  the  porous  anode  of 
SOFC  reconstructed  via  GROMACS,  but  also  thermal  properties 
(including  TC,  VHC,  and  TEC)  and  size  effect  of  the  porous  anode  are 
numerically  calculated  through  Large-scale  Atomic/Molecular 
Massively  Parallel  Simulator  (LAMMPS)  package,  an  acronym  for 
LAMMPS.  In  addition,  the  combined  AA-CG-MD  method  is  devel¬ 
oped  in  detail  and  validated  by  comparing  the  prediction  results 
with  the  open  literature  data,  and  it  is  also  providing  some  insight 
into  an  interconnection  of  sintering  conditions  of  GB  involving  Ni 
and  YSZ  nanocrystals,  the  nanostructures  of  the  porous  anode  in 
morphology  and  the  resulted  thermal  properties  of  the  porous 
anode  at  the  nanoscale. 

2.4.  Numerical  solution  method 

In  this  work,  a  numerical  solution  method  is  developed  to 
reconstruct  and  evaluate  the  porous  anode  based  on  the  above 
three  aspects.  First,  NiO/YSZ  nanocrystals  are  built  at  atomic  level 
using  Atomistic  Simulation  Environment  (ASE)  package,  which 
provides  PYTHON  modules  for  manipulating  atoms  and  deter¬ 
mining  interactions  of  atoms  in  terms  of  DFT.  When  the  atomic 
numbers  and  lattice  constants  are  specified,  the  DFT  calculations  of 
NiO/YSZ  nanocrystals  are  performed  by  the  Jacapo  calculator, 
which  is  an  ASE  interface  for  Dacapo  and  fully  compatible  with  ASE 
package.  In  the  Jacapo,  the  ultra-soft  pseudopotentials  are  used  to 
describe  the  interactions  between  valence  and  core  electron,  and 
the  wave-functions  are  expanded  in  plane  waves  with  the  kinetic 
energy  cutoff  of  350  eV.  To  describe  electron  exchange-correlation 
interactions,  one  adopted  a  generalized  gradient  approximation 
with  the  Perdew-Wang  function  [41  ].  Thus,  the  crystalline  struc¬ 
ture  of  NiO/YSZ  nanocrystals  can  be  created  accurately  without 
information  lost.  In  result,  the  NiO  nanocrystal  is  modeled  by  125 
FCC  lattices,  and  each  NiO  lattice  has  two  nested  FCC  lattices  of 
nickel  and  oxygen.  Similarly,  the  YSZ  nanocrystal  is  also  built  by  125 
tetragonal  lattices.  Secondly,  NiO/YSZ  nanocrystals  are  coarse¬ 
grained  from  atom  representation  to  groups  of  atoms  (beads)  by 
VMD  package  to  reduce  degrees  of  freedom  and  improve  the 
computational  efficiency.  Meanwhile,  interactions  of  the  beads  are 
also  established  using  the  SBCG  module  of  VMD  package.  Neigh¬ 
boring  beads  are  connected  by  harmonic  springs,  while  separate 
molecules  interact  through  nonbonded  forces.  Interactions  are 
parameterized  on  the  basis  of  the  AA  model  and  available  experi¬ 
mental  data.  Therefore,  there  are  around  6  NiO/YSZ  lattices  in  one 
bead  of  the  NiO/YSZ  nanocrystal.  By  this  method,  the  typical  NiO/ 
YSZ  nanocrystal  with  the  size  of  5  x  5  x  5  lattices  is  represented  by 


20  nonpolar  beads,  respectively.  Thirdly,  a  target  box  with  cubic 
50  nm  representing  the  interface  region  between  Ni/NiO  and  YSZ 
phases  is  generated  in  GROMACS  package,  and  then  26,260  NiO 
beads  and  20,000  YSZ  beads  are  randomly  dispersed  into  the  target 
box.  The  energy  minimization  is  performed  in  the  course  of  a  step 
integration  procedure  to  find  a  local  potential  energy  minimum 
near  the  starting  structure.  This  short  energy  minimization  looses 
the  overlapped  beads  at  the  beginning.  MD  running  is  following 
energy  minimization,  which  usually  includes  compressing,  heating 
and  equilibrating  the  system,  respectively.  In  addition,  it  is  often 
advisable  to  perform  equilibration  using  weak-coupling  techniques 
for  temperature  and  pressure,  especially  if  the  system  is  far  from 
equilibration.  The  Berendsen  algorithm  is  applied  here  to  control 
the  system  pressure  and  temperature.  The  equilibrating  process  is 
conducted  subsequently  for  around  200  ns  in  a  canonical  (NTV) 
ensemble.  During  the  simulation,  the  results  and  data  are  auto- 
saved  every  1000  steps  (10  ns)  and  used  for  later  analysis.  In 
result,  the  nanostructure  of  the  anode  is  reconstructed  after  the  MD 
running.  Based  on  the  anodic  nanostructure,  thermal  properties 
can  be  evaluated  systematically  using  LAMMPS  package.  In  this 
study,  all  simulations  are  programmed  in  the  command-line 
interface  from  GROMACS  and  LAMMPS  package.  However,  there 
is  a  format  transformation  from  GROMACS  data  to  an  input  file  of 
LAMMPS,  which  is  implemented  using  MATLAB  code. 

3.  Results  and  discussion 

In  general,  a  solid  oxide  fuel  cell  is  composed  of  four  functional 
layers,  i.e.,  anode,  electrolyte,  cathode  and  interconnect.  The  first 
three  of  these  are  YSZ  ceramic-based.  In  this  work,  a  popular  setup 
of  Ni/YSZ— YSZ— LSM/YSZ  is  employed.  It  should  be  noted  that  the 
ion  conductivity  in  YSZ  is  strongly  dependent  on  the  operating 
temperature.  Consequently,  SOFC  has  to  be  operated  at  appro¬ 
priate  temperatures  ranging  from  800  I<  to  1200  K.  Heat  transfer 
process  has  been  focused  on  in  particular,  which  always  takes 
place  along  with  internal  reforming  reactions  in  the  porous  anode 
for  hydrocarbons,  e.g.,  methane,  ethane,  ethanol,  etc.,  and  elec¬ 
trochemical  reactions,  typically  involving  oxidations  of  H2/CO  in 
the  porous  anode  and  reduction  of  02  in  the  porous  cathode.  On 
the  one  hand,  as  a  feedstock  for  SOFC,  mainly  including  H2,  CO, 
and  CH4,  they  must  be  heated  to  a  certain  temperature,  in  which 
CH4  is  internally  reformed  into  H2/CO  with  heat  energy  consumed 
due  to  the  endothermic  steam  reforming  reaction.  On  the  other 
hand,  reduction  of  oxygen  into  oxygen  ions  takes  place  in  the 
cathode  with  heat  energy  released  due  to  the  exothermic  reaction, 
and  then  the  oxygen  ions  are  transported  through  the  electrolyte 
to  the  anode  where  they  can  electrochemically  oxidize  H2/CO  with 
heat  released  owing  to  the  exothermic  reactions.  In  addition, 
electric  resistance  of  the  materials  also  accounts  for  the  heat 
production  as  a  result  of  Joule  Effect  (JE).  In  other  words,  heat 
production  and  consumption  occurring  in  SOFC,  especially  in  the 
anode,  must  be  managed  dynamically,  because  heat  transfer  can 
either  have  an  effect  on  the  comprehensive  efficiency  of  SOFC,  or 
probably  result  in  thermal  failure.  Therefore,  when  heat  transfer 
process  inside  SOFC  is  mentioned,  some  remarkable  parameters 
must  be  taken  into  account,  as  highlighted  in  the  following 
sessions. 

3.1.  Nanostructure 

To  investigate  effects  on  the  functional  components  of  SOFC 
anodes,  the  sintering  processes  and  nanocrystalline  structures  have 
been  studied  extensively  over  the  past  decades.  In  fabrication 
processes,  the  porous  anode  is  multilayered  co-sintered  in  order  to 
obtain  better  densification  and  dispersity.  Actually,  this 
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manufacturing  process  can  be  simulated  very  well  by  using  the  AA- 
CG-MD  method,  which  totally  mimics  the  crystallization  of  the 
NiO/YSZ  nanoparticles  and  the  self-assembly  phenomena  of  above 
particles  taking  place  in  the  sintering  process.  This  is  quite  inspiring 
for  exploring  the  comprehensive  performance  of  the  porous  anode 
via  the  numerical  method  instead  of  the  expensive  and  complicated 
experiments.  In  this  work,  a  small  domain  (see  Fig.  1)  with  cubic 
50  nm  is  selected  from  all  interface  regions  between  Ni  and  YSZ 
clusters,  which  is  corresponding  to  TPB  features,  as  shown  in  Fig.  2, 
and  then  the  MD  simulation  is  implemented  in  the  thermostat  and 
barostat  conditions,  respectively,  until  the  simulation  system  ach¬ 
ieves  the  equilibration  that  all  components  of  the  porous  anode  can 
be  expressed  with  balanced  porosity,  connectivity,  dispersity  and 
tortuosity.  In  the  simulation  process,  the  Berendsen  weak  coupling 
method  is  applied  in  regulating  the  temperature  and  pressure  of 
the  system  based  on  initially  isotropic  temperature  and  pressure. 
Because  at  the  nanoscale,  atomic  or  molecular  motion  is  dynamic 
and  disordered,  therefore,  it  is  assumed  that  thermal  properties 
relevant  to  the  porous  anode  in  this  study  are  isotropic.  In  addition, 
a  comparison  of  the  MD-based  nanostructure  and  the  TEM-based 
image  of  the  porous  anode  is  discussed  as  follows. 

Fig.  1  shows  numerical  reconstruction  of  the  porous  anode 
obtained  through  the  AA-CG-MD  method  using  GROMACS  pack¬ 
age,  mimicking  the  experimental  fabrication  process.  First,  the  NiO 
and  YSZ  anodic  nanocrystals  are  generated  as  the  sintering  pow¬ 
ders,  including  5x5x5  NiO  lattices  consistent  with  the  exper¬ 
imental  result  as  highlighted  in  Fig.  2  and  5x5x5  YSZ  lattices, 
using  ASE  package,  which  provides  PYTFION-based  modules  for 
building  and  manipulating  atoms,  analyzing  simulations,  visuali¬ 
zations,  etc.,  at  the  atom  level.  Secondly,  the  NiO  and  YSZ  nano¬ 
crystals  are  coarse-grained  into  “pseudo  particles”  (beads), 
representing  groups  of  atoms  by  VMD  package,  which  supplies  a 
SBCG  tool  for  transforming  nanostructures  from  AA  representation 
to  CG  beads.  According  to  a  neural  network  learning  algorithm,  the 
SBCG  approach  is  employed  to  determine  placements  and 


Fig.  1.  The  nanostructure  of  the  porous  anode  of  SOFC  with  cubic  50  nm,  in  which  the 
green  color  represents  NiO  phase,  the  yellow  color  YSZ  phase,  and  the  pores  gas  phase. 
(For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred 
to  the  web  version  of  this  article.) 


Fig.  2.  The  TEM-based  image  of  the  porous  anode  of  SOFC,  in  which  the  white  color 
represents  pores  (gas  phase),  the  silver  color  NiO/Ni  phases,  and  the  gray  color  YSZ 
phase  [20]. 


interactions  of  all  CG  beads,  and  also  to  attach  correct  mass  to  each 
bead.  Thirdly,  26,260  NiO  beads  and  20,000  YSZ  beads  in  total  are 
randomly  dispersed  into  the  50  x  50  x  50  nm3  simulation  box, 
and  then  a  MD  running  is  conducted  over  20  ns  using  GROMACS 
package.  After  the  MD  running,  the  final  structure  configuration  is 
obtained  for  the  equilibrium  state  as  shown  in  Fig.  1.  In  view  of  the 
XY  plane,  an  even  self-assembly  phenomenon  is  exhibited.  The 
NiO  beads  in  green  agglomerate  together  and  compose  the  NiO 
phase  for  conducting  electrons.  The  YSZ  beads  in  yellow  are  con¬ 
nected  to  transport  oxygen  ions.  Pores  are  also  dispersed  evenly 
for  diffusing  gas  and  water.  In  addition,  all  three  phases  are  linked, 
as  expected,  to  provide  the  active  sites  where  the  electrochemical 
reaction  takes  place.  Figs.  1  and  2  show  the  common  morphology 
of  the  porous  anode  of  SOFC,  which  is  basically  dependent  on  the 
sintering  technology.  Meanwhile,  Fig.  1  displays  more  detailed 
nanostructure  of  TPB  regions  in  the  porous  anode  and  the  self- 
assembly  phenomenon  compared  to  Fig.  2  [20].  In  the  begin¬ 
ning,  the  NiO  phase  exists  inside  the  anode,  when  the  electro¬ 
chemical  reactions  start,  the  NiO  phase  will  be  reduced  into  the  Ni 
phase  in  the  hydrogen  atmosphere.  In  fact,  the  Ni  phase  instead  of 
NiO  phase  is  functional  for  conducting  electrons  and  catalyzing  the 
reactions  in  the  operating  process  of  SOFC.  When  the  electro¬ 
chemical  reactions  end,  the  Ni  phase  will  be  oxidized  into  the  NiO 
phase  again.  So  the  Ni  and  NiO  phases  co-exist  in  the  porous 
anode,  and  the  redox  process  between  them  circulates.  In  addi¬ 
tion,  each  phase  must  be  continuous  as  much  as  possible  in  order 
to  conduct  charges  (oxygen-ions  and  electrons)  and  diffuse  gases 
(FI2/CO  and  vapor).  Meanwhile,  all  phases  must  be  dispersed  as 
much  as  possible  because  abundant  TPBs  are  needed  for  electro¬ 
chemical  reactions  taking  place  at  a  high  rate.  According  to  Fig.  1,  it 
is  easily  found  that  the  nano-sintering  technology  contributes  to 
make  the  functional  components  dispersed  very  well,  which  is 
quite  meaningful  for  increasing  the  TPB’s  length.  In  other  words, 
the  self-assembly  phenomena  take  place  at  the  nanoscale, 
inspiringly,  by  which  the  continuity,  dispersity  and  porosity  are 
balanced  very  well  to  some  extent. 
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3.2.  Thermal  conductivity 

Heat  conduction  always  happens  within  or  through  a  body  by 
microscopic  diffusion  and  collisions  of  involved  particles,  e.g., 
electrons,  atoms,  molecules  and  phonons,  etc.,  due  to  a  tempera¬ 
ture  gradient.  In  this  study,  the  porous  anode  is  totally  made  up  of 
NiO  and  YSZ  species  in  the  form  of  CG  particles.  As  a  result,  heat 
conduction  is  mediated  by  the  combination  of  vibrations  and  col¬ 
lisions  of  adjacent  free  electrons,  and  of  propagation,  collisions  and 
scattering  of  phonons.  TC  is  an  inherent  property  of  thermal- 
conductive  materials,  which  is  basically  temperature-dependent 
and  commonly  included  in  Fourier’s  law  for  heat  conduction,  pre¬ 
senting  that  the  local  heat  flux  is  proportional  to  the  negative 
temperature  gradient  as  follows. 


q  =  —kVT  (11) 

So  in  this  case  it  seems  quite  essential  to  estimate  TC  of  the 
anode  of  SOFC.  For  the  engineering  fabrication  of  the  cermet-based 
anode,  a  multi-layered  [42,43]  co-sintering  44,45]  technique  is 
traditionally  employed  to  create  the  porous  and  multi-layered 
anode  from  NiO  and  YSZ  powders,  which  is  always  implemented 
at  high  temperatures  but  below  the  melting  points  of  the  used 
materials.  This  is  so  because  at  higher  temperatures,  phase  trans¬ 
formation  for  YSZ  component  occurs  from  the  monoclinic  crystal 
phase  to  the  tetragonal  crystal  phase.  The  tetragonal  phase  for  YSZ 
nanocrystals  contributes  to  eliminate  the  potential  crack  in  the 
porous  anode,  which  is  prone  to  take  place  due  to  unmatched 
thermal  expansion  between  Ni  and  YSZ  materials.  In  principle, 
appropriately  high  pressures  are  also  applied  at  the  initial  stage  for 
sintering  in  order  to  pursue  the  multi-layered  anode  with  suitable 
porosities  [46,47],  normally  around  30-40  vol%.  In  these  cases,  TC 
is  numerically  calculated  for  the  porous  anodes  reconstructed  at 
various  pressures  ranging  from  5  bar  to  10  bar  and  at  a  certain 
temperature  of  1600  K  according  to  an  equilibrium  MD  approach, 
i.e.,  the  Green— Kubo  (GI<)  method  19,48,49].  In  this  method,  the 
average  TC  in  the  porous  anodes  is  related  to  the  ensemble  average 
of  the  auto-correlation  of  the  heat  flux  into  or  out  of  the  computing 
domain  and  a  universal  temperature. 
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k  =  (/x(°)7x(t)>dt  =  j  (J(0)-J(t))dt  (13) 

0  0 

where  J  is  the  heat  flux,  e*  is  the  per-atom  energy,  including  po¬ 
tential  and  kinetic  energy,  S;  is  the  per-atom  stress  tensor,  V;  is  a 
matrix-vector  multiply. 

As  shown  in  Fig.  3,  TC  of  the  porous  anodes  with  the  composi¬ 
tion  of  Nio.33YSZ  is  slightly  proportional  to  measuring  tempera¬ 
tures,  which  is  accounted  for  by  phonon  heat  conduction. 
Theoretically,  TC  involved  in  phonons  is  usually  described  by  the 
Boltzmann  equation  [34,50,51],  in  which  TC  is  positively  related  to 
measuring  temperature. 
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Fig.  3.  TC  of  the  porous  anodes  with  Ni0.33YSZ  reconstructed  at  various  pressures  and 
1600  K. 
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where  (n)  is  the  phonon  number,  v  is  the  energy  transport  velocity 
of  phonons,  t  is  a  relaxation  time  approximation,  Al  is  the  phonon 
thermal  conductivity,  A  is  the  mean  free  path  for  phonon,  and  O/OTs 
denotes  the  heat  capacity. 

Meanwhile,  it  is  not  difficult  to  notify  that  TC  along  the  vertical 
direction  is  consistent  with  sintering  pressure  for  the  porous  an¬ 
odes.  In  general,  the  sintering  pressure  at  properly  high  levels  can 
easily  make  big  pores  collapse  into  small  pores.  As  a  result,  TC, 
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Fig.  4.  TC  of  the  porous  anodes  reconstructed  at  various  Ni  contents,  10  bar  and 
1600  K. 
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Table  1 

The  comparison  of  predicted  TC  with  open  literature  experi¬ 
mental  data. 

Ref.  [12]  11  W  m-1  K-1 

Ref.  [22]  23  W  m"1  K1 

Ref.  [52]  6.23  W  m"1  K1 

Current  work  2-13  W  m^1  K1 


mechanical  strength  and  TPB’s  length  in  the  high-pressure  sintered 
anodes  are  benefited  simultaneously. 

In  Fig.  4,  the  effect  of  Ni  content  on  TC  in  the  anode  is  taken  into 
account.  The  change  of  TC  in  horizontal  direction  is  similar  to  Fig.  3. 
Vertically,  the  square-symbolized  curve  represents  the  case  with 
the  lowest  Ni  content,  and  it  is  found  that  phonon  heat  conduction 
is  dominating  in  the  corresponding  anode  due  to  relatively  abun¬ 
dant  YSZ  phase.  In  contrast,  the  heat  energy  is  mainly  carried  by 
free  electrons  from  the  Ni  phase  in  both  other  cases,  which  in¬ 
dicates  that  TC  tends  to  flatten  with  increasing  temperature. 

In  both  figures,  the  quantities  of  TC  are  in  good  agreement  with 
the  open  literature  data  [12,22,52],  as  shown  in  Table  1. 

Based  on  Table  1,  the  simulation  results  are  validated  against 
experimental  data,  and  presented  here  to  give  a  detailed  insight 
about  anodic  TC.  It  is  found  that  a  TC  range  from  2  to  13  W  itT1  K_1 
for  the  porous  anodes  is  relatively  reliable,  and  the  difference  can 
be  explained  by  various  compositions  and  sintering  conditions 
applied,  as  discussed  above. 


application  at  a  wide  range  from  room  temperature  to  1200  K.  As  a 
remarkable  advantage  of  sintering,  the  thermal  expansion  of  the 
NiO  and  YSZ  species  can  be  controlled  by  firing  to  generate  crys¬ 
talline  distributions  that  will  influence  the  overall  expansion  of  the 
resulting  material  in  the  desired  direction.  In  addition,  a  compre¬ 
hensive  trade-off  must  be  done  with  an  eye  on  other  properties 
affected,  e.g.,  density,  porosity,  tortuosity,  robustness,  etc.  In  this 
section,  TEC  is  calculated  for  the  porous  anodes,  which  are  sintered 
at  diverse  conditions  as  previously  described. 

Fig.  5  shows  the  average  TEC  of  the  porous  anodes,  which  are 
reconstructed  at  different  sintering  pressures  and  Ni  contents  in 
MD  technique.  Cases  1-6  illustrate  that  the  sintering  pressure  has  a 
slight  effect  on  the  resulting  TEC.  This  means  TEC  of  the  porous 
anodes  with  a  composition  of  Ni0.33YSZ  is  in  a  weakly  positive 
correlation  with  the  sintering  pressure.  In  theory,  thermal  expan¬ 
sion  generally  decreases  with  increasing  bond  energy,  at  the  same 
time,  the  bond  energy  has  a  negative  tendency  with  increasing 
external  pressure.  In  addition,  the  Ni  phase  usually  represents 
higher  thermal  expansion  in  volume,  typically  39E-6  I<-1  [53,54] 
compared  to  the  YSZ  phase  at  10.6E-6  K-1  roughly.  As  a  result, 
thermal  expansion  of  the  anodes  tends  to  be  worse  if  higher  Ni 
contents  are  attempted.  According  to  the  open  literature  data,  a 
combined  value  of  around  12.5E-6  I<-1  [55]  is  tolerated  in  practical 
applications,  which  corresponds  to  the  ones  with  the  low  Ni  con¬ 
tent  (cases  1—6  in  Fig.  5).  However,  it  is  still  poorly  understood  how 
to  match  thermal  expansions  caused  by  Ni  and  YSZ  phases  to  a  good 
balance. 


3.3.  Thermal  expansion 


In  physics,  thermal  expansion  describes  a  tendency  of  materials 
to  change  in  volume  in  response  to  a  temperature  change,  which 
always  takes  place  when  a  substance  is  heated  or  cooled  due  to  the 
particles  moving,  denoted  by  a  temperature-dependent  variable, 
the  so-called  TEC. 


oty 


1 

V 


(19) 


Accurate  estimation  of  thermal  expansion  in  the  cermet-based 
anode  is  a  key  concern  for  a  wide  range  of  reasons,  and  it  is  evi¬ 
denced  that  the  importance  of  thermal  expansion  in  the  porous 
anode  cannot  be  over-emphasized  by  massively  extensive  research 
work.  Otherwise,  the  cermet-based  anode  is  apt  to  crack  because 
YSZ  is  brittle  and  cannot  tolerate  a  sudden  temperature  change. 
Consequently,  the  Ni  phase  and  the  YSZ  phase  have  to  work  well  in 
consort  and  their  thermal  expansion  must  be  matched  in  the 
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Fig.  5.  TEC  of  the  porous  anodes  reconstructed  at  various  pressures,  Ni  contents  and 
1600  K. 


3.4.  Volumetric  heat  capacity 

In  the  process  of  heat  transfer,  VHC  is  often  used  to  measure 
how  much  internal  energy  (totally  including  kinetic  energy  and 
potential  energy)  could  be  stored  in  a  given  volume  of  a  substance 
when  undergoing  a  given  temperature  change  without  phase- 
change. 

r  _  1  AQ 

v  V  AT 

Approximately,  for  a  given  composite,  like  Ni/YSZ  involved  in 
this  work,  the  VHC  is  directly  proportional  to  the  density  and  SHC  of 
these  materials.  At  the  same  time,  there  exists  a  noticeable  positive 
correlation  between  sintering  pressure  and  the  densities  of  the 
sintered  anodes.  In  heat  transfer,  a  higher  value  of  the  VHC  means 
longer  time  for  the  system  to  reach  equilibration.  In  other  words, 
the  big  VHC  (in  proportion  to  the  density  of  the  resulting  material) 
is  desired  because  it  corresponds  to  a  time  period  required  for  a 
change  in  temperature  occurring  in  the  porous  anode. 
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Fig.  6  illustrates  that  the  predicted  VHC  is  in  positive  correlation 
to  sintering  pressure  and  Ni  content.  In  principle,  the  VHC  corre¬ 
sponding  to  a  given  porous  anode  is  subject  to  the  molar  mass, 
volume,  bulk  density  and  the  number  of  involved  atoms,  in  terms  of 
the  Dulong-Petit  law  [56]  as  shown  in  Equations  (21)  and  (22). 

CM  =  3  R  (21) 


Cv 


_  3 Rp  _  3  kBN 
p  M  V 


(22) 


where  C  is  specific  heat  capacity  (J  kg”1  K”1),  M  is  the  molar  mass 
(g  mol-1),  Cv  is  volumetric  heat  capacity  (J  cm”3  K-1),  p  is  bulk 
density  (g  cm-3),  /<b  is  Boltzmann  constant,  N  is  the  total  number  of 
atoms  involved  in  the  system,  and  R  is  the  gas  constant 
(J  K-1  mol-1). 

It  is  clearly  noticed  in  Fig.  6  that,  to  the  cermet  with  a  specific 
composition,  the  increased  VHC  is  due  to  the  system  densification 
with  increasing  sintering  pressure  (cases  2-9).  In  addition,  at  the 
constant  system  volume,  the  VHC  increases  with  increasing  the 
number  of  Ni  atoms,  i.e.,  increasing  Ni  content  (cases  9—12). 

The  predicted  values  of  VHC  in  Fig.  6  are  comparable  to  the  open 
literature  data  [52,57,58]  as  shown  in  Table  2.  In  Ref.  [52],  an 
experimental  value  of  1.49  was  reported  when  the  anode  with  30% 
porosity  was  employed.  A  difference  compared  to  this  work  can  be 
explained  by  the  distinct  anodic  compositions.  Yang  et  al.,  indicated 
a  value  of  1.35  in  their  research  work  [57].  In  Ref.  [58],  the  effect  of 
yttria  content  on  the  VHC  of  YSZ  was  measured  experimentally. 
One  found  that  the  VHC  of  YSZ  decreases  with  increasing  yttria 
content.  In  general,  the  VHC  is  dependent  on  the  anodic  composi¬ 
tion,  as  well  as  Ni  content  and  sintering  pressure.  Unfortunately,  it 
is  a  bit  difficult  to  justify  all  numerical  values  due  to  insufficient 
experimental  data.  However,  the  approaches  applied  in  this  study 
have  been  validated,  and  it  would  be  reasonable  to  mimic  the 
sintering  process  and  predict  relevant  thermal  properties. 


3.5.  Size  effect  on  thermal  conductivity 
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Fig.  7.  Size  effect  on  TC  of  the  porous  anode  rebuilt  at  8  bar  and  1600  K. 


related  to  a  temperature-dependent  parameter,  i.e.,  mean  free  path 
of  phonon  as  shown  in  Equation  (23). 

k  =  1 CmvI  (23) 

where  v  is  the  phononic  velocity  (m  s”1),  l  is  mean  free  path  of 
phonon  (m). 

According  to  another  research  work  [21  ],  it  is  described  that  the 
extrapolation  for  the  bulk  TC  is  relatively  reliable  when  the  system 
size  used  in  the  MD  simulations  is  comparable  to  the  largest  mean 
free  path  of  phonons,  which  can  be  determined  according  to 
Equation  (24). 


kBT 

V2iz  d2p 


(24) 


where  /<b  is  the  Boltzmann  constant  in  J  K_1,  T  is  the  temperature  in 
K,  p  is  pressure  in  Pa,  and  d  is  the  diameter  of  the  involved  particles 
in  m. 


A  lot  of  investigations  for  SOFC  at  the  nanoscale  have  been 
conducted  during  past  years,  which  are  attempted  to  exploit  the 
inspiring  performance  for  SOFC  as  deep  as  possible.  At  the  same 
time,  the  advent  of  the  micro  revolution  also  requires  that  the 
physical  properties  of  involved  materials  must  be  properly  esti¬ 
mated,  which  in  particular  requires  a  suitable  size  of  the  simulation 
domain.  The  effect  of  the  simulation  domain  size  on  the  TC  should 
be  taken  into  account  in  the  thermal  design  for  SOFC  components. 
In  this  study,  the  bulk  material  TC  is  predicted  in  the  porous  anode 
sintered  at  8  bar  and  1600  K  by  using  a  GK  method  in  MD  simu¬ 
lations.  It  is  found  that  a  size-independent  TC  can  be  achieved  with 
the  box  size  no  larger  than  200  nm  for  the  high  temperature  case 
(1198  K  in  Fig.  7). 

In  Fig.  7,  it  is  illustrated  that  predicted  TC  of  the  porous  anode  is 
becoming  varied  with  increasing  size  of  the  simulation  box.  When 
the  box  size  is  larger  than  200  nm,  the  TC  starts  to  increase  more  or 
less.  In  addition,  the  TC  of  the  porous  anode  in  Fig.  7  presents  a 
slight  vibration  relating  to  measuring  temperatures.  Something 
will  account  for  this  phenomenon  that  in  phonon  engineering,  TC  is 

Table  2 

The  comparison  of  predicted  VHC  with  the  open  literature  data. 

Ref.  [52]  1.49  J  cm"3  K"1 

Ref.  [57]  1.35  J  cm"3  K"1 

Ref.  [58]  2.46-3.85  J  cm"3  K1 

Current  work  0.4—5  J  cm-3  K1 


4.  Conclusions 

With  development  in  sciences  and  technologies,  a  rather  big 
amount  of  investigations  for  SOFC  at  the  nanoscale  are  appearing. 
On  the  basis  of  this  background,  the  microscopic  structures  of  the 
porous  anode  of  SOFC  are  reconstructed,  and  then  its  thermal 
properties  are  systematically  calculated  based  on  the  molecular 
dynamics  simulations.  In  this  study,  the  all-atom  coarse-grained 
molecular  dynamics  method  is  developed,  which  contributes  to  the 
above  numerical  computations.  In  general,  when  the  nano¬ 
sintering  technology  is  applied  in  practice  to  pursue  the  good 
performance  of  the  porous  anode  relevant  to  the  well-balance  of 
continuity,  dispersity,  conductivity,  porosity  and  tortuosity,  the  key 
point  is  how  to  control  the  sintering  pressure  and  temperature, 
particle  size  and  composition.  Compared  to  the  conventional 
fabrication  process,  the  nano-sintering  technology  can  produce  the 
well-distributed  nanostructure,  and  a  phenomenon  of  phase  ag¬ 
gregation  is  also  avoided  successfully.  Thermal  conductivity  is  in 
positive  correlation  to  the  sintering  pressure,  and  in  low  Ni  anodes, 
it  slightly  increases  with  measuring  temperature  due  to  phonon 
heat  conduction.  Instead,  in  rich  Ni  anodes,  this  tendency  is  flat¬ 
tening  owing  to  heat  conduction  contributed  by  free  electrons. 
Thermal  expansion  increases  with  increasing  sintering  pressure 
and  Ni  content  due  to  increasing  bond  energy.  The  increased  VHC  is 
due  to  the  system  densification  with  increasing  sintering  pressure. 
In  addition,  at  the  constant  system  volume,  the  VHC  increases  with 
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increasing  Ni  content.  The  effect  of  the  simulation  domain  size  on 
thermal  conductivity  is  also  evaluated.  It  shows  that  the  extrapo¬ 
lation  for  the  bulk  thermal  conductivity  is  reliable  when  the  box 
size  is  comparable  to  the  biggest  mean  free  path  of  phonons 
involved  in  the  simulation  system,  otherwise,  an  oscillation  is 
presented.  Meanwhile,  the  simulation  results  are  compared  against 
experimental  data,  with  reasonable  agreement.  It  is  believed  that  a 
successful  thermal  design  for  the  porous  anode  represents  a 
comprehensive  balance  for  its  sintering  pressure,  temperature,  as 
well  as  particle  size  and  species  composition  in  the  fabrication 
process. 
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